#!/bin/bash

# ABACUS势能面扫描脚本 - Be原子薄膜生长
# 菱形柱状区域扫描，同步执行确保任务顺序完成

# 设置ABACUS路径和运行参数
ABACUS_CMD="/home/zxg/software/abacus3.9.014/abacus-develop/build/abacus"
CORE=10

# 菱形底面顶点坐标（埃）
A=(0 3.21)
B=(2.77994 4.815)
C=(2.77994 1.605)
D=(5.55988 3.21)

# 扫描参数（埃）
Z_START=4.2
Z_END=10.2
STEP=0.4

# 创建结果文件
RESULT_FILE="potential_energy_surface.dat"
echo "# X Y Z Energy Force_X Force_Y Force_Z" > $RESULT_FILE

# 计算扫描范围
min_x=0
max_x=5.55988
min_y=1.605
max_y=4.815

echo "扫描范围: X[$min_x, $max_x], Y[$min_y, $max_y], Z[$Z_START, $Z_END]"
echo "扫描步长: $STEP 埃"

# 计算总点数
total_points=0
for x in $(seq $min_x $STEP $max_x); do
    for y in $(seq $min_y $STEP $max_y); do
        # 检查点是否在菱形内
        # 使用向量叉积方法判断点是否在菱形内
        AB=($(echo "${B[0]} - ${A[0]}" | bc -l) $(echo "${B[1]} - ${A[1]}" | bc -l))
        AP=($(echo "$x - ${A[0]}" | bc -l) $(echo "$y - ${A[1]}" | bc -l))
        BC=($(echo "${C[0]} - ${B[0]}" | bc -l) $(echo "${C[1]} - ${B[1]}" | bc -l))
        BP=($(echo "$x - ${B[0]}" | bc -l) $(echo "$y - ${B[1]}" | bc -l))
        CD=($(echo "${D[0]} - ${C[0]}" | bc -l) $(echo "${D[1]} - ${C[1]}" | bc -l))
        CP=($(echo "$x - ${C[0]}" | bc -l) $(echo "$y - ${C[1]}" | bc -l))
        DA=($(echo "${A[0]} - ${D[0]}" | bc -l) $(echo "${A[1]} - ${D[1]}" | bc -l))
        DP=($(echo "$x - ${D[0]}" | bc -l) $(echo "$y - ${D[1]}" | bc -l))
        
        cross_AB_AP=$(echo "${AB[0]} * ${AP[1]} - ${AB[1]} * ${AP[0]}" | bc -l)
        cross_BC_BP=$(echo "${BC[0]} * ${BP[1]} - ${BC[1]} * ${BP[0]}" | bc -l)
        cross_CD_CP=$(echo "${CD[0]} * ${CP[1]} - ${CD[1]} * ${CP[0]}" | bc -l)
        cross_DA_DP=$(echo "${DA[0]} * ${DP[1]} - ${DA[1]} * ${DP[0]}" | bc -l)
        
        # 如果所有叉积同号，则点在菱形内
        if (( $(echo "($cross_AB_AP >= 0 && $cross_BC_BP >= 0 && $cross_CD_CP >= 0 && $cross_DA_DP >= 0) || ($cross_AB_AP <= 0 && $cross_BC_BP <= 0 && $cross_CD_CP <= 0 && $cross_DA_DP <= 0)" | bc -l) )); then
            for z in $(seq $Z_START $STEP $Z_END); do
                total_points=$((total_points + 1))
            done
        fi
    done
done

echo "总扫描点数: $total_points"
echo "开始扫描..."

# 实际扫描
completed_points=0
for x in $(seq $min_x $STEP $max_x); do
    for y in $(seq $min_y $STEP $max_y); do
        # 检查点是否在菱形内
        AB=($(echo "${B[0]} - ${A[0]}" | bc -l) $(echo "${B[1]} - ${A[1]}" | bc -l))
        AP=($(echo "$x - ${A[0]}" | bc -l) $(echo "$y - ${A[1]}" | bc -l))
        BC=($(echo "${C[0]} - ${B[0]}" | bc -l) $(echo "${C[1]} - ${B[1]}" | bc -l))
        BP=($(echo "$x - ${B[0]}" | bc -l) $(echo "$y - ${B[1]}" | bc -l))
        CD=($(echo "${D[0]} - ${C[0]}" | bc -l) $(echo "${D[1]} - ${C[1]}" | bc -l))
        CP=($(echo "$x - ${C[0]}" | bc -l) $(echo "$y - ${C[1]}" | bc -l))
        DA=($(echo "${A[0]} - ${D[0]}" | bc -l) $(echo "${A[1]} - ${D[1]}" | bc -l))
        DP=($(echo "$x - ${D[0]}" | bc -l) $(echo "$y - ${D[1]}" | bc -l))
        
        cross_AB_AP=$(echo "${AB[0]} * ${AP[1]} - ${AB[1]} * ${AP[0]}" | bc -l)
        cross_BC_BP=$(echo "${BC[0]} * ${BP[1]} - ${BC[1]} * ${BP[0]}" | bc -l)
        cross_CD_CP=$(echo "${CD[0]} * ${CP[1]} - ${CD[1]} * ${CP[0]}" | bc -l)
        cross_DA_DP=$(echo "${DA[0]} * ${DP[1]} - ${DA[1]} * ${DP[0]}" | bc -l)
        
        # 如果所有叉积同号，则点在菱形内
        if (( $(echo "($cross_AB_AP >= 0 && $cross_BC_BP >= 0 && $cross_CD_CP >= 0 && $cross_DA_DP >= 0) || ($cross_AB_AP <= 0 && $cross_BC_BP <= 0 && $cross_CD_CP <= 0 && $cross_DA_DP <= 0)" | bc -l) )); then
            for z in $(seq $Z_START $STEP $Z_END); do
                
                completed_points=$((completed_points + 1))
                progress=$(echo "scale=2; $completed_points * 100 / $total_points" | bc -l)
                echo "进度: $progress% - 扫描点 ($(printf "%.2f" $x), $(printf "%.2f" $y), $(printf "%.2f" $z))"
                
                # 转换为Direct坐标（注意：晶格常数已经包含在晶格矢量中）
                X_dir=$(echo "scale=10; $x / 1.889726 / 11.1198" | bc -l)
                Y_dir=$(echo "scale=10; $y / 1.889726 / 9.63" | bc -l)
                Z_dir=$(echo "scale=10; $z / 1.889726 / 20.0" | bc -l)
                
                # 创建临时STRU文件
                cat > STRU_temp << EOF
ATOMIC_SPECIES
Be    9.012  Be_ONCV_PBE-1.0.upf

NUMERICAL_ORBITAL
Be_gga_10au_100Ry_6s2p.orb

LATTICE_CONSTANT
1.889726

LATTICE_VECTORS
11.11980000000  0.000000000000  0.000000000000
0.000000000000  9.630000000000  0.000000000000
0.000000000000  0.000000000000  20.00000000000

ATOMIC_POSITIONS
Direct

Be
0.000
26
0.250000000000  0.166667000000  0.000000000000  1  1  1  mag  0.0
0.000000000000  0.666667000000  0.000000000000  1  1  1  mag  0.0
0.000000000000  0.000000000000  0.000000000000  1  1  1  mag  0.0
0.250000000000  0.500000000000  0.000000000000  1  1  1  mag  0.0
0.000000000000  0.333333000000  0.000000000000  1  1  1  mag  0.0
0.250000000000  0.833333000000  0.000000000000  1  1  1  mag  0.0
0.333333000000  0.000000000000  0.130325000000  1  1  1  mag  0.0
0.083333000000  0.500000000000  0.130325000000  1  1  1  mag  0.0
0.083333000000  0.166667000000  0.130325000000  1  1  1  mag  0.0
0.333333000000  0.666667000000  0.130325000000  1  1  1  mag  0.0
0.333333000000  0.333333000000  0.130325000000  1  1  1  mag  0.0
0.083333000000  0.833333000000  0.130325000000  1  1  1  mag  0.0
0.750000000000  0.166667000000  0.000000000000  1  1  1  mag  0.0
0.500000000000  0.666667000000  0.000000000000  1  1  1  mag  0.0
0.500000000000  0.000000000000  0.000000000000  1  1  1  mag  0.0
0.750000000000  0.500000000000  0.000000000000  1  1  1  mag  0.0
0.500000000000  0.333333000000  0.000000000000  1  1  1  mag  0.0
0.750000000000  0.833333000000  0.000000000000  1  1  1  mag  0.0
0.833333000000  0.000000000000  0.130325000000  1  1  1  mag  0.0
0.583333000000  0.500000000000  0.130325000000  1  1  1  mag  0.0
0.583333000000  0.166667000000  0.130325000000  1  1  1  mag  0.0
0.833333000000  0.666667000000  0.130325000000  1  1  1  mag  0.0
0.833333000000  0.333333000000  0.130325000000  1  1  1  mag  0.0
0.583333000000  0.833333000000  0.130325000000  1  1  1  mag  0.0
0.000000000000  0.333333000000  0.260650000000  1  1  1  mag  0.0
$(printf "%.10f" $X_dir)  $(printf "%.10f" $Y_dir)  $(printf "%.10f" $Z_dir)  1  1  1  mag  0.0
EOF
                
                # 使用临时STRU文件
                cp STRU_temp STRU
                
                # 同步运行ABACUS计算（等待完成）
                echo "正在计算点 ($(printf "%.2f" $x), $(printf "%.2f" $y), $(printf "%.2f" $z))..."
                OMP_NUM_THREADS=1 mpirun -n $CORE $ABACUS_CMD > out.log 2> err.log
                
                # 检查计算是否成功完成
                if [ $? -eq 0 ] && [ -f "OUT.ABACUS/running_scf.log" ]; then
                    # 提取能量
                    ENERGY=$(grep "FINAL_ETOT_IS" OUT.ABACUS/running_scf.log | tail -1 | awk '{print $2}')
                    
                    # 提取最后一个原子的力（新添加的原子）
                    if [ -f "OUT.ABACUS/FORCE_IS.txt" ]; then
                        # 获取最后一行力数据（第26个原子）
                        FORCE_LINE=$(tail -n 1 OUT.ABACUS/FORCE_IS.txt)
                        FORCE_X=$(echo $FORCE_LINE | awk '{print $2}')
                        FORCE_Y=$(echo $FORCE_LINE | awk '{print $3}') 
                        FORCE_Z=$(echo $FORCE_LINE | awk '{print $4}')
                    else
                        FORCE_X="N/A"
                        FORCE_Y="N/A"
                        FORCE_Z="N/A"
                    fi
                    
                    # 写入结果文件
                    echo "$x $y $z $ENERGY $FORCE_X $FORCE_Y $FORCE_Z" >> $RESULT_FILE
                    echo "计算完成: 能量 = $ENERGY, 力 = ($FORCE_X, $FORCE_Y, $FORCE_Z)"
                else
                    echo "$x $y $z ERROR ERROR ERROR ERROR ERROR" >> $RESULT_FILE
                    echo "警告: 点 ($(printf "%.2f" $x), $(printf "%.2f" $y), $(printf "%.2f" $z)) 计算失败"
                fi
                
                # 清理临时文件，保留结果
                rm -f STRU_temp out.log err.log
                rm -rf OUT.ABACUS
                
            done
        fi
    done
done

echo "扫描完成！结果保存在 $RESULT_FILE"
echo "总共完成 $completed_points/$total_points 个点的计算"